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optimization, focusing on the convergence of the average run as a function of run- 
time and system size. The method has a single free parameter, which we determine 
numerically and justify using a simple argument. Our numerical results demonstrate 
that on random graphs, extremal optimization maintains consistent accuracy for in- 
creasing system sizes, with an approximation error decreasing over runtime roughly 
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the average run suggests that these are far from optimal, with large fluctuations 
between individual trials. But when only the best runs are considered, results con- 
sistent with theoretical arguments are recovered. 
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I. INTRODUCTION 

Optimizing a system of many variables with respect to some cost function is a task 
frequently encountered in physics. The determination of ground-state configurations in 
disordered materials [1, 2, 3, 4] and of fast-folding protein conformations [5] are but two 
examples. In cases where the relation between individual components of the system is frus- 
trated [6], the cost function often exhibits a complex "landscape" [7] in configuration space, 
posing challenges to neighborhood search procedures. Indeed, for growing system size the 
cost function may exhibit a rapidly increasing number of unrelated local extrema, separated 
by sizable barriers that can make the search for the exact optimal solution unreasonably 
costly. It is of great importance to develop fast and reliable approximation methods for 
finding optimal or acceptable near-optimal solutions with high probability. 

In recent papers we have introduced a new method, called extremal optimization (EO), 
to tackle such hard optimization problems [8, 9]. EO is a method based on the dynamics of 
non-equilibrium processes and in particular those exhibiting self-organized criticality [10], 
where better solutions emerge dynamically without the need for parameter tuning. Previ- 
ously, we have discussed the basic EO algorithm, its origin, and its performance compared 
with other methods. We have demonstrated that the algorithm can be adapted to a wide 
variety of NP-hard problems [11]. We have shown that for the graph partitioning problem, a 
simple implementation of EO yields state-of-the-art solutions, even for systems of iV > 10^ 
variables [8]. For large graphs of low connectivity, EO has been shown to be faster than 
genetic algorithms [12] and more accurate than simulated annealing [13], two other widely 
applied methods. A numerical study [14] has shown that EO's performance relative to sim- 
ulated annealing is particularly strong in the neighborhood of phase transitions, "where the 
really hard problems are" [15]. In fact, preliminary studies of the phase transition in the 
3-coloring problem [16] as well as studies of ground state configurations in spin glasses [3, 17] 
suggest that EO may become a useful tool in the exploration of low-temperature properties 
of disordered systems. 

In the present work we focus on the intrinsic features of the method, by investigating its 
average performance. For this purpose, we have performed an extensive numerical study 
of EO on the graph bipartitioning problem. We have considered various kinds of graph 
ensembles, both with geometric and with random structure, for an increasing number of 
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vertices N. The results show that for random graphs, EO converges towards the optimal 
configuration in a power-law manner, typically requiring no more than 0{N) update steps. 
For geometric graphs the averaged large- A?" results are less convincing, but if we instead focus 
on the best out of several trials, near-optimal results emerge. Our implementation of EO 
has one single tunable parameter, and we find a simple relation to estimate that parameter 
given the allowed runtime and system size. Many of our numerical results here have been 
independently confirmed by J. Dall [18]. 

The paper is organized as follows. In Sec. II we introduce the graph bipartitioning 
problem, and in Sec. Ill we describe the extremal optimization algorithm. Sec. IV deals in 
detail with our numerical results. In Sec. V we conclude with an outlook on future work. 

II. GRAPH BIPARTITIONING 
A. Definition 

The graph bipartitioning problem (GBP) is easy to formulate. Take N vertices, where 
N is an even number, and where certain of the vertex pairs are connected by an edge. Then 
divide the vertices into two sets of equal measure N/2 such that the number of edges con- 
necting both sets, the "cutsize" m, is minimized. The global constraint of an equal division 
of vertices places the GBP among the hardest problems in combinatorial optimization, since 
determining the exact solution with certainty would in general require a computational efi^ort 
growing faster than any power of [19]. It is thus important to find "heuristic" methods 
that can obtain good approximate solutions in polynomial time. Typical examples of ap- 
phcations of graph partitioning are the design of integrated circuits (VLSI) [20] and the 
partitioning of sparse matrices [21]. 

The general description of a graph in the previous paragraph is usually cast in more spe- 
cific terms, defining an ensemble of graphs with certain characteristic. These characteristics 
can affect the optimization problem drastically, and often refiect real-world desiderata such 
as the geometric lay-out of circuits or the random interconnections in matrices. Therefore, 
let us consider a variety of different graph ensembles, some random and some geometric in 
structure. 
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B. Classes of graphs studied 

One class of graphs that has been studied extensively is that of random graphs without 
geometric structure [22]. Here, edges between any two vertices are taken to exist with 
probability p: on the average, an instance has a total of = pN{N — l)/2 vertices and 
the mean connectivity per vertex is a = p{N — 1). Following standard terminology we refer 
to graphs of this sort as the ensemble of random graphs, even though the other classes of 
graphs we consider all have stochastic properties as well. 

Another often-studied class of graphs without geometric structure is generated by fixing 
the number of connections a at each vertex, but with random connections between these 
vertices [23, 24]. In particular, we consider here the case where a — 3: the ensemble of 
trivalent graphs, randomly connected graphs with exactly three edges originating from each 
vertex. 

The third class we consider is an ensemble with geometric structure, where the vertices 
are situated on a cubic lattice. Edges are placed so as to connect some (but not all) nearest 
neighbors on the lattice: a ratio x of the total number of nearest-neighboring pairs are 
occupied by an edge, and those edges are distributed at random over the possible pairs. 
For a cubic lattice, the average connectivity is then given hy a = 6x. This class of graphs 
corresponds to a dilute ferromagnet, where each lattice site holds a ±-spin and some (but not 
all) nearest-neighboring spins possess a coupling of unit strength. Here, the GBP amounts 
to the equal partitioning of -|- and — spins while minimizing the interface between the two 
types [25], or simply finding the ground state under fixed (zero) magnetization. We thus 
refer to this class as the ensemble of ferromagnetic graphs. 

The final class we consider is that of geometric graphs specified by N randomly distributed 
vertices in the 2-dimensional unit square, where we place edges between all pairs whose two 
vertices are separated by a distance of no more than d [26]. The average connectivity is then 
given by q; = Nird^. The GBP on this class of graphs has the advantage of a simple visual 
representation, shown in Fig. 1. Again following standard terminology, we refer to this class 
simply as the ensemble of geometric graphs. 

It is known that graphs without geometric structure, such as those in the first two classes, 
are typically easier to optimize than those with geometric structure, such as those in the final 
two classes [26]. The characteristics of the GBP for non-geometric and geometric graphs 
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FIG. 1: Plot of a geometric graph with N = 500 vertices and average connectivity a = 6, par- 
titioned into 250 red and 250 blue vertices. Starting from an initial random assignment of red 
and blue colors, EO arrives at near-optimal configurations consisting of domains of red and blue 
vertices, separated by an interface across which "bad" edges (green lines) connect poorly-adapted 
vertices. 

at low connectivity appear to be very different, due to the dominance of long loops in the 
former and short loops in the latter. The ensemble of random graphs has a structure that is 
locally tree-like, allowing for a mean- field treatment that yields some exact results [25]. By 
contrast, the ensemble of geometric graphs corresponds to continuum percolation of "soft" 
(overlapping) circles, for which precise numerical results exist [27]. 

Each of the graph ensembles that we consider is characterized by a control parameter, 
the average connectivity a. The difficulty of the optimization problem for each type varies 
significantly with a. In this study we focus on sparse graphs for which a is kept constant, 
independent of A^. Sparse graphs have very different properties from the dense graphs studied 
by Fu and Anderson [28]. These sparse graphs are generally considered to pose the most 
difficult partitioning problems, and extremal optimization (EO) is particularly competitive 
in this regime [14] . In order to facilitate our average performance study, we fix a to a given 



6 



value on each ensemble. For random graphs, where the connectivity varies between vertices 
according to a Poisson distribution, letQ;=p(A^ — 1)=2 (and so p ~ 1/-^)- For trivalent 
graphs, by construction a = 3. For ferromagnetic graphs, let a = = 2. For geometric 
graphs, let a = NudP = 6. In all of these cases, the connectivity is chosen to be just above 
the phase transition at ctcrit, below which the cutsize m almost always vanishes [14]. These 
critical regions are especially interesting because they have been conjectured to coincide 
with the hardest-to-solve instances in many combinatorial optimization problems [15, 29]. 

Finally, in light of the numerous comparisons in the physics literature between the GBP 
and the problem of finding ground states of spin glasses [1], it is important to point out 
the main difference. This is highlighted by the ensemble of ferromagnetic graphs. Since 
couplings between spins are purely ferromagnetic, all connected spins invariably would like 
to be in the same state; there is no possibility of local frustration. Frustration in the GBP 
merely arises from the global constraint of an equal partition, forcing spins along an interface 
to attain an unfavorable state (see Fig. 1). All other spins reside in bulk regions where they 
can maintain the same state as their neighbors. In a spin glass, on the other hand, couplings 
can be both ferromagnetic and anti-ferromagnetic. Spins everywhere have to compromise 
according to conflicting conditions imposed by their neighbors; frustration is local rather 
than global. 

C. Basic scaling arguments 

If we neglect the fact that the structure of these sparse graphs is that of percolation 
clusters, we can obtain some elementary insights into the expected scaling behavior of the 
optimal cutsize with increasing size N, m ^ N^/'^ . For graphs without geometric structure 
(random graph ensemble and trivalent graph ensemble), one can expect that the cutsize 
should grow linearly in i.e., v — 1. Indeed, this argument can be made rigorous for 
arbitrary connectivity a. Extremal optimization performs very well on these graphs, and 
previous numerical studies using EO all give i/ f» 1 [14] . 

For graphs with geometric structure (ferromagnetic graph ensemble and geometric graph 
ensemble), the value of v is less clear. We can approximate a graph with a ci-dimensional 
geometric structure as a hyper-cubic lattice of length L = N'^^'^, where the lattice sites are 
the vertices of the graph and the nearest-neighbor bonds are the edges, of which only a finite 
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fraction are occupied. There are thus about E N edges in the graph. To partition the 
graph, we are roughly looking for a (d — 1) -dimensional hyper-plane cutting the graph into 
two equal-sized sets of vertices. Such an interface between the partitions would cut ~ L'^~^ 
bonds, and thus ~ N^~^/'^ edges. Following this argument, the 3-d ferromagnetic graphs 
should have a cutsize scaling with N'^^^ and the 2-d geometric graphs should have a cutsize 
scahng with N^^"^. 

However, while this may the case for a typical partition of the graph, it may not be 
the case for an optimal partition. The interface for an optimal cut of a sparse graph could 
well be much rougher than our argument suggests, taking advantage of large voids between 
clusters of connected vertices. The number of cut edges would then be much below the 
estimate based on assuming a flat interface, making l/v < 1 — 1/d. In our previous studies 
using EO, we found 1/u ^ 0.75 ± 0.05 for ferromagnetic graphs and l/z/fti0.6±0.1for 
geometric graphs [14], i.e., above the upper bound, and our newer results do not improve 
on these (seen later in Fig. 6). This could indicate that the actual values are close to the 
upper bound, but also that for graphs with geometric structure EO fails to find the optima 
on instances of increasing size. 

Similar behavior has been observed with other local search methods [26], reflecting the fact 
that sparse geometric graphs generally pose a much greater challenge than do sparse random 
graphs. In contrast, a heuristic such as METIS [30], a hierarchical decomposition scheme 
for partitioning problems, works better for geometric graphs than for random graphs [31]. 
METIS performs particularly well for sparse geometric graphs, and typically produces better 
results than EO for a = 6. Furthermore, if speed is the dominant requirement, METIS is 
superior to any local search method by at least a factor of N. But for random graphs at 
q; = 2 or the trivalent graphs, METIS' results are poor compared to EO's, and for all type 
of graphs METIS' performance deteriorates with increasing connectivity. 

III. EXTREMAL OPTIMIZATION ALGORITHM 
A. Motivation 

The extremal optimization method originates from insights into the dynamics of non- 
equilibrium critical phenomena. In particular, it is modeled after the Bak-Sneppen mecha- 
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nism [32], which was introduced to describe the dynamics of co-evolving species. 

Species in the Bak-Sneppen model are located on the sites of a lattice, and each one has 
a "fitness" represented by a value between and 1. At each update step, the smallest value 
(representing the most poorly adapted species) is discarded and replaced by a new value 
drawn randomly from a flat distribution on [0, 1]. Without any interactions, all the fitnesses 
in the system would eventually approach 1. But obvious interdependencies between species 
provide constraints for balancing the system's overall condition with that of its members: the 
change in fitness of one species impacts the fitness of an interrelated species. Therefore, at 
each update step, the Bak-Sneppen model replaces the fitness values on the sites neighboring 
the smallest value with new random numbers as well. No explicit definition is provided for 
the mechanism by which these neighboring species are related. Yet after a certain number 
of updates, the system organizes itself into a highly correlated state known as self-organized 
criticality (SOC) [10]. In that state, almost all species have reached a fitness above a certain 
threshold. But these species merely possess what is called punctuated equilibrium [33] : since 
only one's weakened neighbor can undermine one's own fitness, long periods of "stasis", 
with a fitness above the threshold, are inevitably punctuated by bursts of activity. This co- 
evolutionary activity cascades in a chain reaction ( "avalanche" ) through the system. These 
fiuctuations can involve any number of species, up to the system size, making any possible 
configuration accessible. Due to the extremal nature of the update, however, the system as 
a whole will always return to states in which practically all species are above threshold. 

In the Bak-Sneppen model, the high degree of adaptation of most species is obtained by 
the elimination of poorly adapted ones rather than by a particular "engineering" of better 
ones. While such dynamics might not lead to as optimal a solution as could be engineered 
under specific circumstances, it provides near-optimal solutions with a high degree of latency 
for a rapid adaptation response to changes in the resources that drive the system. A sim- 
ilar mechanism, based on the Bak-Sneppen model, has recently been proposed to describe 
adaptive learning in the brain [34] . 

B. Algorithm description 

Inspired by the Bak-Sneppen mechanism, we have devised the EO algorithm with the goal 
of accessing near-optimal configurations for hard optimization problems using a minimum 
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of external control. Previously, we have demonstrated that the EO algorithm is apphcable 
to a wide variety of problems [11, 17]. Here, we focus on its implementation for the GBP. 
In the GBP, EO [8] considers each vertex of a graph as an individual variable with its 

own fitness parameter. It assigns to each vertex i a "fitness" 



where gi is the number of "good" edges connecting i to other vertices within its same set, 
and bi is the number of and "bad" edges connecting i to vertices across the partition. (For 
unconnected vertices we fix Aj = 1.) Note that vertex i has an individual connectivity of 
Oil — Qi + bi, while the overall mean connectivity of a graph is given by a = ZliCti/-^ and 
the outsize of a configuration is given hj m = J2ibi/2. 

At all times an ordered list is maintained, in the form of a permutation 11 of the vertex 
labels i such that 



and i — Ii{k) is the label of the kth ranked vertex in the hst. 

Feasible configurations have N/2 vertices in one set and N/2 in the other. To define a 
local search of the configuration space, we must define a "neighborhood" for each configu- 
ration within this space [35]. The simplest such neighborhood for the GBP is given by an 
exchange of (any) two spins between the sets. With this exchange at each update we can 
search the configuration space by moving from the current configuration to a neighboring 
one. In close analogy with the Bak-Sneppen mechanism, our original parameter-free imple- 
mentation of EO merely swapped the vertex with the worst fitness, 11(1), with a random 
vertex from the opposite set [8]. Over the course of a run with tmax update steps, the cutsize 
of the configurations explored varies widely, since each update can result in better or worse 
fitnesses. Proceeding as with the gap-equation in the Bak-Sneppen model [36] , we can define 
a function ■m{t) to be to be the cutsize of the best configuration seen during this run up to 
time t. By construction m{t) is monotonically decreasing, and m{t„iax.) is the output of a 
single run of the EO- algorithm. 

We find that somewhat improved results are obtained with the following one-parameter 
implementation of EO. Draw two integers, I < ki,k2 < N, from a probability distribution 



^9i/{9i + bi), 



(1) 



An(i) < An(2) < • • • < An(Ar), 



(2) 



P{k) (X k~\ {l<k< N), 



(3) 
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on each update, for some r. Then pick the vertices ii — Il{ki) and ^2 = n(/c2) from the 
rank-ordered hst of fitnesses in Eq. (2). (We repeatedly draw k2 until we obtain a vertex in 
the opposite set from ki.) Let vertices ii and i2 exchange sets no matter what the resulting 
new outsize may be. Then, reevaluate the fitnesses A for ii, i2, and all vertices they are 
connected to {2a on average). Finally, reorder the ranked list of A's according to Eq. (2), 
and start the process over again. Repeat this procedure for a number of update steps per run 
that is linear in system size, imax = AN, and store the best result generated along the way. 
Note that no scales to limit fluctuations are introduced into the process, since the selection 
follows the scale-free power-law distribution P{k) in Eq. (3) and since — unlike in heat bath 
methods — all moves are accepted. Instead of a global cost function, the rank-ordered list of 
fitnesses provides the information about optimal configurations. This information emerges 
in a self-organized manner merely by selecting with a bias against badly adapted vertices, 
rather than ever "breeding" better ones. 



C. Discussion 

1. Definition of fitness 

We now discuss some of the finer details of the algorithm. First of all, we stress that we use 
the term "fitness" in the sense of the Bak-Sneppen model, in marked contrast to its meaning 
in genetic algorithms. Genetic algorithms consider a population of configurations and assign 
a fitness value to an entire configuration. EO works with only a single configuration and 
makes local updates to individual variables within that configuration. Thus, it is important 
to reiterate that EO assigns a fitness Aj to each of the system's N variables, rather than to 
the system as a whole. 

While the definition of fitness in Eq. (1) for the graph partitioning problem seems natural, 
it is by no means unique. In fact, in general the sum of the fitnesses does not even represent 
the cost function we set out to optimize, because each fitness is locally normalized by the 
total number of edges touching that vertex. It may seem more appropriate to define fitness 
instead as Xi — gi, the number of "good" connections at a vertex, or else as Aj = —bi, which 
amounts to penalizing a vertex for its number of "bad" connections. In both cases, the 
sum of all the fitnesses is indeed linearly related to the actual cost function. The first of 
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these choices leads to terrible results, since almost all vertices in near-optimal configurations 
have only good edges, and so in most cases gi is simply equal to the connectivity of the 
vertex. The second choice does lead to a viable optimization procedure and one that is 
easily generalizable to other problems, as we have shown elsewhere [11]. But for the GBP, 
we find that the results from that fitness definition are of poorer quality than those we 
present in this paper. It appears productive to consider all vertices in the GBP on an equal 
footing by normalizing their fitnesses by their connectivity as in Eq. (1), so that Aj e [0, 1]. 
The fact that each vertex's pursuit towards a better fitness simultaneously minimizes its 
own contribution to the total cost function ensures that EO always returns sufficiently close 
to actual minima of the cost function. 

Note that ambiguities similar to that of the fitness definition also occur for other optimiza- 
tion methods. In general, there is usually a large variety of different neighborhoods to choose 
from. Furthermore, to facihtate a local neighborhood search, cost functions often have to 
be amended to contain penalty terms. It has long been known [26] that simulated annealing 
for the GBP only becomes effective when one allows the balanced partition constraint to be 
violated, using an extra term in the cost function to represent unbalanced partitions in the 
cost function. Controlling this penalty term requires an additional parameter and additional 
tuning, which EO avoids. 

2. The parameter r 

Indeed, there is only one parameter, the exponent r in the probability distribution in 
Eq. (3), governing the update process and consequently the performance of EO. It is intuitive 
that a value of r should exist that optimizes EO's average-case performance. If r is too small, 
vertices would be picked purely at random with no gradient towards any good partitions. If 
T is too large, only a small number of vertices with particularly bad fitness would be chosen 
over and over again, confining the system to a poor local optimum. Fortunately, we can 
derive an asymptotic relation that estimates a suitable value for r as a function of the allowed 
runtime and the system size N. The argument is actually independent of the optimization 
problem under consideration, and is based merely on the probability distribution in Eq. (3) 
and the ergodicity properties that arise from it. 

We have observed numerically that the choice of an optimal r coincides with the transition 
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of the EO algorithm from ergodic to non-ergodic behavior. But what do we mean by 
"ergodic" behavior, when we are not in an equihbrium situation? Consider the rank-ordered 
hst of fitnesses in Eq. (2), from which we choose the individual variables to be updated. If 
r = 0, we choose variables at random and can reach every possible configuration of the 
system with equal probability. EO's behavior is then perfectly ergodic. Conversely, if r 
is very large, there will be at least a few variables that may never be touched in a finite 
runtime t, because they are already sufficiently fit and high in rank k. Hence, if there 
are configurations of the system that can only be reached by altering these variables first, 
EO will never explore them and accordingly is non-ergodic. Of course, for any finite r, 
different configurations will be explored by EO with difi^crcnt probabilities. But we argue 
that phenomenologically, we may describe EO's behavior as ergodic provided every variable, 
and hence every rank on the list, gets selected at least once during a single run. There will 
be a value of r at which certain ranks, and therefore certain variables, will no longer be 
selected with finite probability during a given runtime. We find that this value of r at the 
transition to non-ergodic behavior corresponds roughly with the value at which EO displays 
its best performance. Clearly, this makes the choice of r dependent on the runtime and on 
the system size A^. 

Assuming that this coincidence indicates a causal relation between the ergodic transition 
and optimal performance, we can estimate the optimal r in terms of runtime t and size N. 
EO uses a runtime tmax = -^N, where A is a constant, probably much larger than 1 but much 
smaller than N. (For a justification of the runtime scahng hnearly in N, see Sec. JVC 2.) 
We argue that we are at the "edge of ergodicity" when during AN update steps we have a 
chance of selecting even the highest rank in EO's fitness list, k = N, about once: 

P{k = N)AN ~ 1. (4) 

With the choice of the power-law distribution for P{k) in Eq. (3), we obtain 

{t -1)N-^AN r^l, (N^oo), (5) 

where the factor of r — 1 arises from the norm of P{k). Asymptotically, we find 

log (^/ log AT) (g) 

logA^ ^ ^ ^ ^ 

Of course, large higher-order corrections may exist, and there may well be deviations in the 
optimal T among different classes of graphs since this argument does not take into account 
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graph structure or even the problem at hand. Nevertheless, Eq. (6) gives a qualitative 
understanding of how to choose r, indicating for instance that it varies very slowly with N 
and A but will most likely be significantly larger than its asymptotic value of unity Not 
surprisingly, with the numerical values A ^ 10^ and N ^ 10'^ used in previous studies, we 
typically have observed optimal performance for r ~ 1.3 — 1.6 (see also Rcf. [18]). Our 
numerical study of r is discussed in Sec. IV B below. (We note that in this study we often 
use runtimes with A > N to probe the extreme long-time convergence behavior of EO. In 
that case, Eq. (6) can not be expected to apply. Yet, the optimal value of r still increases 
with A, as will be seen in the numerical results in Sec. IV B and Fig. 3a.) 

3. Efficient ranking of the fitness values 

Strictly speaking, the EO-algorithm as wc have described it has a cost proportional to 
2aN^logN per run. One factor of N arises simply from the fact that the runtime, i.e., 
the number of update steps per run, is taken to scale linearly with the system size. The 
remaining factor of 2Q;A^logA^ arises from the necessity to maintain the ordered list of 
fitnesses in Eq. (2): during each update, on average 2a vertices change their fitnesses and 
need to be reordered, since the two vertices chosen to swap partitions arc each connected on 
average to a other vertices. The cost of sequentially ordering fitness values is in principle 
NlogN. However, to save a factor of N, we have instead resorted to an imperfect heap 
ordering of the fitness values, as already described in Ref. [8]. Ordering a hst of N numbers 
in a binary tree or "heap" ensures that the smallest fitness will be at the root of the tree, 
but does not provide a perfect ordering between all of the other members of the list as a 
sequential ordering would. Yet, with high probability smaller fitnesses still reside in levels 
of the tree closer to the root, while larger fitnesses reside closer to the end-nodes of the tree. 
To maintain such a tree only requires the 0{logN) moves needed to replace changing fitness 
values. 

Specifically, consider a hst of N fitness values. This hst will fill a binary tree with at most 
^max+1 levels, where Zmax — [log2(-^)] ([^] denotes the integer part of x) and I — 0,1, ... , /max, 
where / = is the level consisting solely of the root, / = 1 is the level consisting of the two 
elements extending from the root, etc. In general, the Ith level contains up to 2^ elements, 
and all level are completely filled except for the end-node-level Imax, which will only be 
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partially filled in case that N < 2'™^''+^ — 1. Clearly, by definition, every fitness on the Zth 
level is worse than its two descendents on the (l + l)th level, but there could nevertheless 
be fitnesses on the Ith level that are better than some of the other fitnesses on the (/ + l)th 
(or even greater) level. On average, though, the fitnesses on the Ith level are always worse 
than those on the (/ + l)th level, and better than those on the (/ — l)th level. 

Thus, instead of applying the probability distribution in Eq. (3) directly to a sequentially 
ordered list, we save an entire factor of N in the computational cost by using an analogous 
(exponential) probability distribution on the (logarithmic) number of levels I in our binary 
tree, 

Q{1) (X 2-(--i)', (0 < Z < [log,{N)] + 1). (7) 

Prom that level I we then choose one of the 2' fitnesses at random and update the vertex 
associated with that fitness value. Despite the differences in the implementations, our studies 
show that heap ordering and sequential ordering produce quite similar results. In particular, 
the optimal value of r found for both methods is virtually indistinguishable. But with the 
update procedure using the heap, EO algorithm runs at a cost of merely 0(2Q;iVlog A^). 

Under some circumstances, it may be possible to maintain a partially or even perfectly 
ordered list at constant cost, by using a hash table. For instance, for trivalent graphs or 
more generally for ct-valent graphs, each vertex i can be in only one of a+1 attainable states, 
bi — 0, 1, . . . , a and Aj = bi/a. Thus, instead of time consuming comparisons between A's, 
fitness values can be hashed into and retrieved from "buckets" , each containing all vertices 
with a given b^. For an update, we then obtain ranks according to Eq. (3), determine which 
bucket that rank points to, and retrieve one vertex at random from that bucket. 

Even in cases where the fitness values do not fall neatly into a discrete set of states, such a 
hash table may be an effective approximation. But great care must be taken with respect to 
the distribution of the A's. This distribution could look dramatically different in an average 
configuration and in a near-optimal configuration, because in the latter case fitness values 
may be densely clustered about 1. 
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4- Startup routines 

The results obtained during a run of a local search algorithm can often be refined using 
an appropriate startup routine. This is an issue of practical importance for any optimization 
method [35]. For instance, Ref. [26] has explored improvements for the partitioning of geo- 
metric graphs by initially dividing the vertices of the graph into two geometrically defined 
regions (for instance, drawing a line through the unit square) . This significantly boosted the 
performance of the Kernighan-Lin algorithm on the graphs. Such methods are not guaran- 
teed to help, however: simulated annealing shows little improvement using a startup [26]. 
Happily, the performance of the EO algorithm typically is improved considerably with a 
clever startup routine. 

Previously, we have explored a startup routine using a simple clustering algorithm [8], 
which accomplishes the separation of the graph into domains not only for geometric but also 
for random graphs. The routine picks a random vertex on the graph as a seed for a growing 
cluster, and recursively incorporates the boundary vertices of the cluster until N/2 vertices 
are covered, or until no boundary vertices exist anymore (signaling a disconnected cluster in 
the graph). The procedure then continues with a new seed among the remaining vertices. 
Such a routine can substantially enhance EO's convergence especially for geometrically de- 
fined graphs (see Sec. IV C 3). In this paper, though, we did not use any startup procedure 
because we prefer to focus on the intrinsic features of the EO algorithm itself. Instead, all 
the results presented here refer to runs starting from random initial partitions, except for a 
small comparison in Sec. JVC 3. 

IV. NUMERICAL RESULTS 

A. Description of EO runs 

In our numerical simulations, we considered the four classes of graphs introduced in 
Sec. II. For each class, wc have studied the performance of EO as a function of the size of 
the graph A^, the runtime t, and the parameter r. To obtain a statistically meaningful test 
of the performance of EO, large values of N were chosen; EO performed too well on smaller 
graphs. The maximum value of N varied with each kind of graph, mostly due to the fact 
that some types required averaging over many more instances, thus limiting the attainable 
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sizes. 

The precise instance sizes were: N — 1022, 2046, 4094 and 8190 for both random and 
trivalent graphs; N = 13^(= 2197), 16^{= 4096) and 20^{= 8000) for ferromagnetic graphs; 
and = 510, 1022, and 2046 for geometric graphs. For each class of graphs, we generated 
a sizable number of instances at each value of N: 32 for random and trivalent graphs, 16 
for ferromagnetic graphs, and 64 for geometric graphs. For each instance, we conducted a 
certain number of EO runs: 8 for random and trivalent graphs, 16 for ferromagnetic graphs, 
and 64 for geometric graphs. Finally, the runtime (number of update steps) used for each 
run was tmax — AN, with A = 512 for random graphs, A = 4096 for trivalent graphs, 
A = 1000 for ferromagnetic graphs and A = 2048 for geometric graphs. 

B. Choosing r 

In previous studies, we had chosen r = 1.4 as the optimal value for all graphs. In light of 
our discussion in Sec. Ill C 2 on how to estimate r, here we have performed repeated runs 
over a range of values of r on the instances above, using identical initial conditions. The goal 
is to investigate numerically the optimal value for r, and the dependence of EO's results on 
the value chosen. In Figs. 2a-d we show how the average cutsizc depends on r, given fixed 
runtime. There is a remarkable similarity in the qualitative performance of EO as a function 
of T for all types and sizes of graphs. Despite statistical fluctuations, it is clear that there is 
a distinct minimum for each class, and that as expected, the results get increasingly worse 
in all cases for smaller as well as larger values of r. 

While the optimal values for r are similar for all types of graphs, there is a definite drift 
towards smaller values of r for increasing A^. Studies on spin glasses [17] have led to similar 
observations, supporting the argument for the scaling of r that we have given in Sec. IIIC2. 
Our data here do not cover a large enough range in N to analyze in detail the dependence 
of r on logN. However, we can at least demonstrate that the results for trivalent graphs, 
where statistical errors are relatively small, are consistent with Eq. (6) above. For fixed N 
but increasing values oi A — t/N, we see in Fig. 3a that the optimal value of r appears 
to increase linearly as a function of log A in the regime 1 <^ A <^ N. At the same time, 
fixing A and increasing A^, we see in Fig. 3b that the optimal value of r appears to decrease 
linearly as a function of 1/logA^, towards a value near t — 1. Thus, for 1 <C -C A^, r 
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(a) Random Graphs at A=51 2 (b) Trivalent Graphs at A=51 2 

_ , ^ , ^ , ^ , ^ r, 0.130 n ' ' ' ' ' ' — 




FIG. 2: Outsize found by EO (in units of N) as a function of r at a fixed value of ^ = t/N ~ 500, 
averaged over all runs on all graphs of a given size N. Results are shown for (a) random graphs, 
(b) trivalent graphs, (c) ferromagnetic graphs, and (d) geometric graphs. For each type of graph 
the minimum shifts very slowly to smaller values of r for increasing A^. 

seems to be converging very slowly towards 1, with a finite correction of ~ log A/ log A^. 
This is in accordance with our estimate, Eq. (6), discussed earlier in Sec. IIIC2. (Note that 
at A ^ 500 in Figs. 2a-b, the optimal values for random and trivalent graphs are close to 
r ^ 1.3, consistent with Eq. (6). In contrast, in our long-time studies below, A N and a 
value of r = 1.45 seems preferable, consistent with Fig. 3a.) 

The foregoing data, including the results plotted in Figs. 2a-d, arise from averages {m{t)) 
over all runs (denoted by (. . .)) and all instances (denoted by an overbar). But it is important 
to note that the conclusions drawn with respect to the optimal value of r from these plots 
are valid only if there is little difference between the average run {'m{t)) and the best run 
'^bcst for a particular instance. While this is the case for the random and trivalent graphs, 
there is a significant difference between {m{t)) and mbest for instances of ferromagnetic and 
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(a) (b) 




5 10 0.05 0.1 

FIG. 3: Plot of Topt for trivalent graphs as a function of (a) log2(^) = 'log2{t/N) for various fixed 
values of N and (b) l/log2(iV) for various fixed values of ^ = t/N. These data points were 
determined by performing runs as in Fig. 2b, and finding the minimum of a quartic fit to data 
obtained at r = 1.2, 1.25, 1.3, . . . , 1.95. In (a) the data increase roughly linearly with log2(A) in 
the regime 1 <^ A <^ N at fixed A^, while in (b) the data extrapolate roughly toward r = 1 for 
^ oo at fixed A, both in accordance with Eq. 6. For ^ ^ 1 this scaling appears to break down, 
while for A N the linear scaling happens to remain valid. Note that the values of Topt in (b) 
correspond to the data points in (a) for log2(^) = 5 and 8. 

geometric graphs, as is discussed later in Sec. JVC 3. In fact, for geometric graphs (Fig. 6 
below), average and best cutsizes often differ by a factor of 2 or 3! Fig. 2d indicates that 
the optimal value of r for the average performance on geometric graphs of size = 2046 is 
below r = 1.2. If we instead plot the fraction of runs that have come, say, within 20% of 
the best value mtest found by EO (which most likely is still not optimal) for each instance, 
the optimal choice for r shifts to larger values, r ^ 1.3 (see Fig. 4). 

We may interpret this discrepancy as follows. At lower values of r, EO is more likely 
to explore the basin of many local minima but also has a smaller chance of descending to 
the "bottom" of any basin. At larger values of r, all but a few runs get stuck in the basin 
of a poor local minimum (biasing the average) but the lucky ones have a greater chance of 
finding the better states. For geometric graphs, the basins seem particularly hard to escape 
from, except at the very low values of r where EO is unlikely to move toward the basin's 
minimum. Thus, in such cases we find that we get better average performance {m{t)) at a 
lower value of r, but the best result mbest of multiple runs is obtained at a higher value of 
r. 
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FIG. 4: Fraction of EO runs on geometric graphs that have come within 20% of the best ever found 
(for each instance) as a function of r, for each value of N. Maxima indicate optimal choice of r for 
finding a few good results among many runs. These maxima occur at higher values of r than the 
minima corresponding to best average performance in Fig. 2. As seen below in Fig. 6, it becomes 
increasingly hard to come close to these best values. 

Clearly, in order to obtain optimal performance of EO within a given runtime and for a 
given class of graphs, further study of the best choice of r would be justified. But for an 
analysis of the scaling properties of EO with runtime and A^, the specific details of how the 
optimal T varies are not significant. Therefore, to simply the following scaling discussion, 
we will simply fix r to a near-optimal value on each type of graph. 

C. Scaling behavior 

1. General results 

As explained in Sec. Ill, the cutsize of the current configuration will fluctuate wildly at 
all times during an EO run and will not in itself converge. Instead, we have to keep track 
of the best configuration obtained so far during a run. Thus, even for times t < tmax during 
an EO run, we refer to the cutsize m{t) of the current best configuration as the "result" of 
that run at time t. Examples of the stepwise manner in which m(t) converges towards the 
optimal cutsize, for a particular instance of each type of graph, are given in Fig. 5. Up until 
times t N, m{t) drops rapidly because each update rectifies poor arrangements created by 
the random initial conditions. For times t ^ A^, it takes collective rearrangements involving 
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(a) Random Graph (b) Trivalent Graph 
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(c) Ferromagnetic Graph (d) Geometric Graph 




FIG. 5: Log-log plot of the outsize as a function of the number of update steps t for a single (large) 
instance of a (a) random, (b) trivalent, (c) ferromagnetic and (d) geometric graph. The solid black 
line represents m{t) for a single run, and the solid red line with error bars represents the average 
{m{t)) over all runs on that instance. In (a) and (b), error bars are very small, indicating that 
there are only small run-to-run fluctuations about {m{t)). By contrast, in (c) and (d), these run- 
to-run fluctuations are huge. In fact, for the geometric graph in (d) we have plotted two separate 
solid black lines representing two selected runs, one poor and one good one. This shows the large 
variations between runs that lead to the large error bars. 

many vertices to (slowly) obtain further improvements. 

While during each run m{t) decreases discontinuously, for the random and trivalent graphs 
the jumps deviate relatively little from that of the mean behavior obtained by averaging over 
many runs ((. . .)). Fluctuations, shown by the error bars in Figs. 5a-b, thus are small and will 
be neglected henceforth. For the ferromagnetic and geometric graphs, these fluctuations can 
be enormous. In Fig. 6 we compare the average performance (m) with the best performance 
mbest foi' each type of graph, at the maximal runtime tmax and averaged over all instances at 
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FIG. 6: Log-log plot of the cutsize (m), averaged over all runs and instances (filled symbol), and 
the best cutsize mbest from all runs on an instance, averaged over all instances (open symbol). This 
is shown for random (o), trivalent (□), ferromagnetic (o), and geometric graphs (A), as a function 
of size A^. The error bars refer only to run-to-run (and not instance-to-instance) fluctuations. For 
random and trivalent graphs, both average and best cutsizes scale linearly in A^, as expected. For 
ferromagnetic and geometric graphs, the best results are several standard deviations better than 
the average results, with a widening gap for increasing N on the geometric graphs. The scaling of 
^best ~ N^/'^ gives u = 1.3 for ferromagnetic graphs, and is consistent with z/ = 2 for geometric 
graphs. 

each A^. The results demonstrate that for random and trivalent graphs the average and best 
results are very close, whereas for ferromagnetic and geometric graphs they are far apart 
(with even the scaling becoming increasingly different for the geometric case). Therefore, 
in the following we will consider scaling of the average results for the first two classes of 
graphs, but properties of the best results for the latter two. 



2. Random and trivalent graphs 

Averaging m(t) for any t over all runs ((...)) and over all instances (overbar), we obtain 
the averaged cutsize 



{m) =m{N,t,T) (8) 

as a function of size N, runtime t, and parameter r. In Fig. 7 we plot m{N, t, r) for 
random and trivalent graphs as a function of t, for each and at a fixed value r = 1.45 
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(a) Random Graphs (b) Trivalent Graph 




FIG. 7: Log-log plot of the average cutsize m{N, t, r) as a function of runtime t at fixed r = 1.45, for 
(a) random graphs and (b) trivalent graphs. The average is taken over runs as well as over instances 
(32 instances for random and 8 instances for trivalent graphs). Error bars represent instance-to- 
instance fluctuations only, and increase with N much more slowly than the mean result. In each 
case, N increases from bottom to top. 

that is near the optimal value for the maximal runtimes tmax (see Sec. IV B). The error 
bars shown here are due only to instance-to-instance fluctuations, and are consistent with 
a normal distribution around the mean. (Note that the error bars are distorted due to the 
logarithmic abscissa.) The fact that the relative errors decreases with N demonstrate that 
self- averaging holds and that we need only focus on the mean to obtain information about 
the typical scaling behavior of an EO run. 

We wish to study the scaling properties of the function m in Eq. (8). First of all, we find 
that generally 

m(iV, t, r) ~ N^'^mit/N, r) (t > > 1), (9) 

reflecting the fact that for EO, as well as for most other graph partitioning heuristics that are 
based on local search, runtimes scale linearly with A^. (This is justified by the fact that each 
of the variables only has 2 states to explore. In, say, the traveling salesperson problem, by 
contrast, each of cities can be reconnected to 0{N) other cities, and so runtimes typically 
scale at least with A^^ [37].) In Fig. 8 we plot m{N,t,T)/N for fixed r = 1.45 as a function 
t/N. We find indeed that the data points from Fig. 7 collapse onto a single scaling curve 
m{t/N, r), justifying the scaling ansatz in Eq. (9) for u = 1. The scaling collapse is a strong 
indication that EO converges in 0{N) updates towards the optimal result, and also that the 
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FIG. 8: Scaling collapse of the data from Fig. 7 onto a single scaling curve m{t/N,T) = m/N as a 
function of t/N at fixed r = 1.45, for (a) random graphs and (b) trivalent graphs. For t > N, data 
points are fitted to the power-law expression in Eq. (10) with numerical values given in Tab. I. 

optimal cutsize itself scales linearly in A^. 

The scaling function m appears to converge at large times to the average optimal result 
(iTT'opt) = {fnopt)/N according to a power law: 

t 



mit/N, t) ~ (mopt) + C (^^) (t » iV » 1). 



(10) 



Fitting the data in Fig. 8 to Eq. (10) for t/N > 1, we obtain for each type of graph a 
sequence of values for (rriopt) and 7 for increasing A^, given in Tab. I. In both cases we find 
values for {rhopt)/N that are quite stable, while the values for 7 slowly decrease with A^. 
The variation in 7 as a function of may be related to the fact that for a fixed r, EO's 



TABLE I: Sequence of values of the fit of Eq. (10) to the data in Fig. 8 for t/N > 1, for each N. 



Graph Type 


Size N 


("T-opt) 


C 


7 


Random 


1022 


0.0423 


0.028 


0.48 




2046 


0.0411 


0.029 


0.46 




4094 


0.0414 


0.032 


0.45 




8190 


0.0425 


0.034 


0.44 


Trivalent 


1022 


0.1177 


0.052 


0.43 




2046 


0.1159 


0.057 


0.41 




4094 


0.1159 


0.062 


0.41 




8190 


0.1158 


0.066 


0.40 



24 



performance at fixed A — t/N deteriorates logarithmically with N, as seen in Sec. IV B. 
Even with this variation, however, the values of 7 for both types of graph are remarkably 
similar: 7 0.4. This implies that in general, on graphs without geometric structure, we can 
halve the approximation error of EO by increasing runtime by a factor of 5-6. This power- 
law convergence of EO is a marked improvement over the mere logarithmic convergence 
conjectured for simulated annealing [13] on NP-hard problems [2]. 

Finally, the asymptotic value of (mopt) obtained for the trivalent graphs can be compared 
with previous simulations [14, 24]. There, the "energy" £ — —1 -|-4(mopt)/3 was calculated 
using the best results obtained for a set of trivalent graphs. Ref. [24] using simulated 
annealing obtained £ = —0.840, while in a previous study with EO [14] we obtained £ = 
—0.844(1). Our current, somewhat more careful extrapolation yields (mopt) = 0.1158A^ 
or £^ = —0.845(1). Even though this extrapolation is based on the average data rather 
than only the best of all runs, the very low fluctuations between runs (see Figs. 5b and 6) 
indicate that the result for S would not change significantly. Thus, the replica symmetric 
solution proposed in Refs. [23, 25] for this version of the GBP, which gives a value of 
£: = -2 X 0.7378/^ = 0.852, seems to be excluded. 

3. Ferromagnetic and geometric graphs 

Unlike on the preceding graphs, EO gives significantly different results for an average run 
and for the best run on geometrically structured graphs (see Figs. 5 and 6). In Fig. 6, at 
least the results from the best run (averaged over all instances) comes close to the scahng 
behavior expected from the considerations in Sec. II C: a fit gives v ^ 1.3 for ferromagnetic, 
and 1/^2 for geometric graphs, while the theory predicts u = 3/2 and u = 2 respectively. 
But even these best cutsizes themselves vary significantly from one instance to another. 
Thus, it is certainly not useful to study the "average" run. Instead, we will consider the 
result of each run at the maximal runtime tmax, extract the best out of k runs, and study 
these results as a function of increasing k. 

Fig. 9 shows the difficulty of finding good runs with near-optimal results for increasing 
size A^. While for ferromagnetic graphs it is possible that acceptable results may be obtained 
without increasing k dramatically, for geometric graphs an ever larger number of runs seems 
to be needed at large N. It does not appear possible to obtain consistent good near-optimal 
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(a) Ferromagnetic graphs (b) Geometric Graphs 
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FIG. 9: Extrapolation plot for the best-of-/c trials for (a) ferromagnetic graphs at r = 1.4 and (b) 
for geometric graphs at r = 1.3. The data for this plot are extracted from the results at tmaxj 
averaging the best-of-Zc results over 16 instances for ferromagnetic, and 64 instances for geometric 
graphs. For comparison, the leftmost data point for each at A: = 16 for ferromagnetic graphs 
and at /c = 64 for geometric graphs corresponds to the "best" results plotted in Fig. 6 for those 
graphs. 

results at a fixed k for increasing A^. 

We saw in Sec. IV C 2 that for random graphs, computational time is well spent on a 
few, long EO runs per instance. We can not address in detail the question of whether, 
for geometrically defined graphs, computational time is better spent on k independent EO- 
runs with t^ax update steps or, say, on a single EO-run with k x t^ax update steps. While 
experience with our data would indicate the former to be favorable, an answer to this 
question depends significantly on and, of course, on the choice of r (see Sec. IV B). Here 
we consider this question merely for a single value, r = 1.3, for which we have run EO 
on the same 64 geometric graphs up to 16 times longer than tmax, but using only k = 4 
restarts. In each of these 4 runs on an instance we recorded the best result seen at multiples 
^ X ^max with n = 1, 2, 4, 8, and 16. For example, the best-of-4 runs at = 1 of this runtime 
corresponds to the best-of-fc results in Fig. 9 for k = 4, while n = 16 would correspond to 
the same amount of running time as k = 64. Fig. 10 shows that fewer but longer runs are 
slightly more successful for larger A^. 

Finally, we have also used the clustering algorithm described in Sec. IIIC4 and Ref. [8] 
on this set of 64 graphs with r = 1.3. For comparison, we again use the best-of-4 runs with 
averages taken at times n x tmax, n = 1,2, and 4. Results at short runtimes improve by a 
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FIG. 10: Equivalent-runtime comparison between different strategies to improve EO results on 
geometric graphs for N = 1022 and 2046 at r = 1.3. The horizontal axis is proportional to the 
inverse of the number of updates used, t = k x n x tmax- Filled symbols refer to the n = 1 results 
for geometric graphs already plotted in Fig. 9b, where k is varied. Open symbols on the dotted 
lines refer to the k = 4 (best-of-4) results, where n is varied. Opaque symbols on the dashed line 
refer to /c = 4 results as well, but using initial conditions generated by a clustering algorithm (see 
Sec. IIIC4). At sufficiently long runtime all strategies are about equal, though with fewer but 
longer runs having a slight edge over more and shorter runs for large N. Even the advantages of 
a non-random initial configuration become less significant at longer runtimes. 

huge amount with such a procedure, but its advantage is eventually lost at longer runtimes. 
While this procedure is cheap, easy, and apparently always successful for geometric graphs, 
our experiments indicate that its effect may be less significant for random graphs and may 
actually result in diminished performance when used for trivalent graphs in place of random 
initial conditions. Clearly, clustering is tailored more toward geometric graphs satisfying a 
triangular inequality. 

V. CONCLUSIONS 

Using the classic combinatorial optimization problem of bipartitioning graphs as an appli- 
cation, we have demonstrated a variety of properties of the extremal optimization algorithm. 
We have shown that for random graphs, EO efficiently approximates the optimal solution 
even at large size N, with an average approximation error decreasing over runtime t as 
t~^'^. For sparse, geometrically defined graphs, finding the ideal (sub-dimensional) interface 




27 



partitioning the graph becomes ever more difficult as size N increases. EO, hke other local 
search methods [26], gets stuck ever more often in poor local minima. However, when we 
consider the best out of multiple runs with the EO algorithm, we recover results close to 
those predicted by theoretical arguments. 

We believe that many of our findings here, notably with regard to EO's fitness definition 
and the update procedure using the parameter r, are generic for the algorithm. Our results 
for optimizing 3-coloring and spin glasses appear to bear out such generalizations [17]. In 
view of this observation, a firmer theoretical understanding of our numerical findings would 
be greatly desirable. The nature of EO's performance derives from its large fluctuations; 
the price we pay for this is the loss of detailed balance, which is the theoretical basis for 
other physically- inspired heuristics, such as simulated annealing [13]. On the other hand, 
unlike in simulated annealing, we have the advantage of dealing with a Markov chain that 
is independent of time [38]. This suggests that our method may indeed be amenable to 
theoretical analysis. We believe that the results EO delivers, as a simple and novel approach 
to optimization, justify the need for further analysis. 
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